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We describe the energy relaxation process produced by surface damping on lattices of classical 
anharmonic oscillators. Spontaneous emergence of localised vibrations dramatically slows down 
dissipation and gives rise to quasi-stationary states where energy is trapped in the form of a gas of 
weakly interacting discrete breathers. In one dimension (ID), strong enough on-site coupling may 
yield stretched-exponential relaxation which is reminiscent of glassy dynamics. We illustrate the 
mechanism generating localised structures and discuss the crucial role of the boundary conditions. 
For two-dimensional (2D) lattices, the existence of a gap in the breather spectrum causes the 
localisation process to become activated. A statistical analysis of the resulting quasi-stationary 
state through the distribution of breathers' energies yield information on their effective interactions. 
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Discrete breathers have been widely analysed 
as mathematical objects, but the relevance of 
their role in physical systems is still debated. In 
fact, any mechanism yielding the spontaneous for- 
mation of breathers is expected to provide inter- 
esting indications in this direction. In this re- 
spect, the phenomenon of relaxation to breather 
states by cooling lattices from their boundaries is 
one of the most interesting examples. Here we 
provide a comprehensive description of this phe- 
nomenon in low dimensional lattices. In partic- 
ular, we review former results and add new in- 
formation, mainly concerning the 2D case. This 
study is based on a combination of numerics and 
theoretical arguments, while we analyse both dy- 
namical and statistical features of the problem to 
clarify the whole scenario. 



I. INTRODUCTION 

Nonlinearity has revealed one of the key ingredients for 
describing many relevant features of different states of 
matter. In the realm of lattice dynamics, scattering pro- 
cesses among phonons, propagation of solitary waves and 
slow energy relaxation are typical examples. Recently, 
considerable efforts have been devoted to the study of 
periodic, localised, non-linear lattice excitations, named 
"breathers" 0] . They are quite peculiar (but generic) 
objects emerging from the interplay of nonlinearity and 
space discreteness Their mathematical properties, 
like existence, stability, mobility etc. have been progres- 
sively unveiled, while they have been identified in many 



different scenarios. 

Their role in non equilibrium dynamics seems to be 
particularly fascinating. An example is the relaxation to 
energy equipartition of short-wavelength fluctuations pj . 
Another interesting scenario where breathers are found 
to emerge spontaneously is observed upon cooling the 
lattice at its boundaries |[ ||, i.e. by considering a 
non-equilibrium process in which energy exchange with 
the environment is much faster at the surface than in the 
bulk. The numerical simulations neatly show that the en- 
ergy dissipation rate may be significantly affected by the 
spontaneous excitation of breathers. As the latter exhibit 
a very weak interaction among themselves and with the 
boundaries, the energy release undergoes a sudden slow- 
ing down and is hardly detected on the time scales of a 
typical simulation. Thus, the lattice remains frozen in 
a pseudo-stationary, metastable configuration which is 
far from thermal equilibrium, which we shall refer to as 
residual state. Despite the absence of disorder, there is 
a close similarity with the glassy behaviour observed in 
disordered systems, as was already pointed out || ||. 

But what are the mechanisms leading to spontaneous 
localisation? As we have shown in a previous paper ||, 
the latter is intimately related to how dissipation acts 
on vibrational modes of different wavelengths. If long- 
wavelength phonons can be efficiently damped out, the 
modulation instability of short lattice waves Q becomes 
highly favoured and breathers can easily emerge from an 
interacting soliton-gas j|, ||]. In this respect, this type 
of non-equilibrium condition is much more effective in 
exciting localised modes than an equilibrium one (see e.g. 
for a related discussion) . 

Spontaneous localisation upon cooling has been ob- 
served for both on-site J|, 0] and pure nearest-neighbour 
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(n.n.) interactions j|, 10 . The two classes of models are 
known to behave differently, e.g. as to the mobility of 
breathers, which is much higher in the absence of local 
coupling. This is also confirmed by relaxations experi- 
ments, where differences in the energy decay laws have 
been observed. In particular, the approach to the resid- 
ual state follows a simple exponential law in the Fermi- 
Pasta-Ulam (FPU) system @, but turns to stretched ex- 
ponential in the presence of local coupling Q . 

Our aim in this paper is to provide an up-to-date 
overview of the phenomenon of breather localisation by 
cooling. Generalities about the energy release process, 
including models and indicators, will be presented in Sec- 
tion [[]], along with the discussion of the cooling process 
in an harmonic lattice and the results obtained for ID 
non linear systems. In particular we shall survey the dif- 
ferent cooling pathways appearing in pure n.n. and on- 
site potentials. A specific subsection will be devoted to 
clarify the interpretation of some recent results |h], [ll| . 
In Section III, we present the results for the n.n. 2D 
FPU model. We want to point out that the 2D is not 
a straightforward extension of the ID case. Indeed, the 
existence of an energy activation threshold for breather 
solutions [Q leads to conjecture that a thermalised state 
may be cooled to the residual state only above some ini- 
tial energy/temperature. Moreover, the residual state is 
a static multi-breather state, whereby in the ID case it 
usually contains a single breather. Finally, we perform 
an analysis of the distribution of breather energies, based 
on a simple statistical model. 



II. RELAXATION AND LOCALISATION IN ID 

In this section we will describe the lattice models and 
the generalities of a typical relaxation experiment as well 
as the relevant indicators used throughout the paper. Let 
us consider a chain of N atoms of unit mass and denote 
with Up the displacement of the p-th particle from its 
equilibrium position pa (a is the lattice spacing). The 
atoms are labelled by the index p = 0, 1, ... N — 1 and 
their dynamics is given by the equations of motion 



= V'(u : 
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where V(x) and U (x) are the interparticle and on-site po- 
tentials, respectively. We assume that V'(0) = U'(0) = 
0, denoting with a prime the derivative with respect to 
x. For convenience, we shall adopt non dimensional units 
such that a and V"(0) are set to unity. Moreover, as we 
want to deal with systems in a finite volume, we will im- 
pose either free-end (w_i = wo, "iv-i = u n) or fixed-end 
(u-i =upf = 0) boundary conditions (BC). 

The last term in eq. ([!]) represents the interaction of 
the atoms with a "zero-temperature" heat bath in the 
form of a linear damping with rate 7. Following j|, |), 
we restricted to the case in which dissipation selectively 



acts only on the atoms located at the chain edges. As 
mentioned above, this choice is crucial and should model 
a physical situation in which energy exchange with the 
environment is much faster at the surface than in the 
bulk. 

The general layout of a simulation can be summarised 
as follows . First, an equilibrium microstate is generated 
by, say, Nose-Hoover (canonical) method pj or by letting 
the Hamiltonian (microcanonical) system (7 = 0) evolve 
for a sufficiently long transient raj. In the following we 
will follow the second strategy. The initial condition for 
the transient is assigned by setting all displacements u p 
to zero and by drawing velocities at random from a Gaus- 
sian distribution. The velocities are then rescaled by a 
suitable factor to fix the desired value of the energy den- 
sity (energy per particle) eo- The resulting set of u p and 
iip is then used as initial condition to integrate Eqs. (|l|) 
with 7 > (see again Rcf. || for details). 

The result of the relaxation process is the residual 
state, i.e. a state with a finite fraction of the initial en- 
ergy stored in the form of breathers. Remarkably, such 
residual state is characterised by a decay time scale which 
is several orders of magnitude longer than that of the lo- 
calisation process. 



A. The harmonic lattice 

A number of features which characterise the relaxation 
dynamics of non-linear lattices can be better understood 
by first studying the harmonic lattice, namely 

V(x) = l -x\ 

We anticipate that the presence of an harmonic on-site 
potential does not alter the main conclusions. We thus 
take [7 = throughout this section. 

For small damping, an approximate analytical solution 
can be found by time-dependent perturbation theory. 
Let 0l = 4sin 2 (<? Q /2) and i] a (a = 0, 1, ... N - 1) denote 
the eigenvalues and normalized eigenvectors of the unper- 
turbed Hamiltonian problem, where the allowed wave- 
numbers are q a = an/N and q a — (a + l)ir/(N +1) for 
free-end and fixed-end BC, respectively. To the lowest 
order in 7 one obtains || 



JV-l 



+ i0a) t 



r\p- 



(1) 



where T] p is the p-th component of the eigenvector r\ a 
and c a are constant amplitudes fixed by the initial con- 
ditions. Notice that, since the number of the damped 
particles is fixed, the perturbative approximation is ex- 
pected to improve upon increasing the system size N. 
The wavenumber-dependent damping rates are found to 



3 



be 

1 J i cos 2 ($f) for free-end BC 
Ta |^sin 2 (<j Q ) for fixed-end BC . 

with To — A/27. The free-end and fixed-end BC sys- 
tems thus show considerably different behaviours. In 
the former case, the least damped modes are the short- 
wavelength ones (a » N), the largest lifetime being 
t jv _ 1 « 2iV 3 /7r 2 7, while the most damped modes are 
the ones in the vicinity of a = 0, with tq being the 



shortest decay time. On the contrary, for fixed ends the 
most damped modes are those around the band center 
(a w A/2) while short- and long-wavelength ones dissi- 
pate very weakly, being t n _ 1 w N(N + l) 2 /2n 2 ^. As we 
shall illustrate in the following, such difference is crucial 
for the localisation in the non-linear system. 

Using the solution ([I]), we evaluated the chain energy 
in the case of equipartition i.e. by replacing c 2 2 with 
its average value 2cq. Finally, recalling eqs. (|2|) and ap- 
proximating for large N the sum over a with an integral 
(t q — » r(q)), we obtained Q 



E(0) 



-2t/r(q) 



dq 



-t/TQ T 
J 



-t/T 
1 



for t <C T , 
for i > r . 



(3) 



where 7 is the modified zero-order Bessel function, 
whose asymptotic expansions have been used to obtain 
the last equality. We conclude that relaxation is exponen- 
tial only at short times while a crossover to a power-law 
decay occurs at t w t . Moreover, despite the differences 
in the relaxation times @ , the decay law for E turns out 
to be the same in the free and fixed-end systems. 

The above formulas must be taken with caution when 
dealing with a large but finite system. In such a case, 
there is of course a lower cutoff at wavenumbers of order 
ir/N. Therefore, after a time of the order of the lifetime 
of the longest-lived mode, t n _ ± ~ N 3 /^, formula (||) no 
longer holds and a further crossover to the exponential 
decay law exp(— 2i/r JV _i) occurs. This has been also ver- 
ified numerically H. 



B. Non— linear lattices: the role of discreteness 

Since the first results of relaxation experiments in non- 
linear lattices were reported [Q ||], there has been some 
debate as to the nature of the decay law of the system 
energy and its relation with the spontaneously emerging 
localised modes. In particular, it has been claimed that 
the phenomenon of energy pinning in the form of discrete 
breathers would result in a glassy-like relaxation, i.e. 
stretched exponential decay law. Here we shall describe 
how the main role in determining the relaxation proper- 
ties of a non-linear lattice is indeed played by the type of 
potential energy In particular, we will highlight the im- 
portance of the relative strength of the inter-particle and 
on-site potentials, i.e. the degree of discreteness of the 
system. In order to do so, we study the general class of 
non-linear lattices obeying the equations of motion (Q). 
In particular, we shall report here the results of numerical 



simulations performed with 

V(x) = -x 2 + -x 4 U{x) = -kx 4 . (4) 
2 4 4 

In this case, the relative strength of local and inter- 
particle non-linearities is accounted for by the parameter 

K. 

In Fig. |l| we show the space-time contour plots of the 
symmetrised site energies, defined as 

h p = ^ii 2 + - [V(u p+ i - u p ) + V{u p - Up_i)] + U(u p ). 

The instantaneous total energy of the system is then 
given by E — J2p=o hp- We see that a single localised 
excitation emerges from the relaxation process for both 
small and large values of n. The pathway to localisation 
is the same as the one observed in the FPU model || . In 
particular, the basic mechanism leading to localisation 
is modulational instability of short wavelength modes. 
The latter can only be effective if dissipation of long 
wavelength phonons is fast enough. As it shows from 
Eq. (||), this occurs only in the case of free-ends BC, 
whereas fixed-end BC strongly inhibit such process. On 
the other hand, breather mobility is strongly reduced in 
the highly discrete system (Fig. |l| (b)). Such difference 
is even more evident in the decay law of the normalized 
energy E(t)/ E(0), which turns from pure exponential to 
stretched-exponential behaviour upon increasing k. This 
is better appreciated by plotting the indicator 

V{t) = -\n[E{t)/E{Q)] (5) 

in a log-log scale, so that a stretched-exponential law of 
the form E(t) = E(0) exp[— (£/r) <T ] becomes a straight 
line with slope a that intercepts the j/-axis at — alnr 
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FIG. 1: Relaxation in an FPU chain with quartic on-site potential, 7 = 0.1, N = 100, eo = 1. Space-time contour plots of the 
symmetrised site energies, (a) k = 0.1. (b) k = 10 
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FIG. 2: FPU chain with quartic on-site potential, N = 100, 7 = 0.1 and eo = 1. (a) Plot of T> vs time. The dashed line is 
the exponential exp(— t/r n ), while the solid line is a fit with the stretched exponential law exp[— (t/r)"] (a ~ 0.61). The arrow 
marks the time at which the system with k — 10 departs from the initial exponential trend, (b) Localisation parameter vs time 
for different values of the on-site coupling k. The data of both T> and C are averaged over 32 initial conditions. 



(Fig. g (a)). The scaling regions correspond to the on- 
set of energy localisation, as can be seen by plotting the 
localisation parameter 



C(t) =N- 



£M0 



(6) 



(Fig. @ (b)). From its very definition, the fewer sites the 



energy is localised onto, the closer C is to N . On the 
other hand, the more evenly the energy is spread among 
all the particles, the closer £ is to a constant of order 
1. It should be noted, that £ is a relative quantity, and 
bears no information on the amount of energy that is 
localised. As a matter of fact, the fraction of eo that gets 
trapped in the residual state is found to be greater the 
higher the degree of discreteness of the system. 

To summarise, if the on-site force is weak enough, the 
behavior is basically the same as the FPU chain which, 
in turn, exhibits in the approach to the residual state the 
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same exponential decay law of its linearised counterpart, 
i.e. E(t)/E{0) oc exp(-i/ T o) (see Eqs. @). The reasons 
for such behaviour are the same discussed for the FPU 
model (g). First of all, an effective harmonic Hamilto- 
nian with energy-dependent renormalized frequencies is 
known to account for several equilibrium properties of 
the FPU chain [jl2| |l3). Second, the harmonic approx- 
imation becomes increasingly accurate as time elapses, 
simply because more and more energy is extracted from 
the system by the reservoir. On the contrary, if the sys- 
tem is highly discrete, the first relaxation stage is indeed 
described by a stretched exponential law. A quantitative 
explanation of it is still lacking, but it is clear that in this 
case an effective description of such genuine non-linear 
phenomenon in terms of quasi-harmonic modes breaks 
down. In other words, the resulting interactions among 
linear modes inhibit the energy flow from the bulk to the 
boundaries, thus slowing down the dissipation. 

The onset of a pseudo-stationary state corresponds in 
Fig. H to the almost flat final portions of the curves. Ac- 
tually, being the system globally dissipative, energy is 
still at that stage exponentially decreasing but with a 
huge time constant 77,. The latter is expected to be in- 
versely proportional to the breather amplitude at the lat- 
tice edges, which, in turn, should be of order exp(— N/t), 
with I <C TV being the localisation length. 

The above scenario is generic for large enough systems. 
Nonetheless, failure to spontaneously localise energy oc- 
curs the more frequently the smaller is the lattice. This 
is because the probability of occurrence of large enough 
energy fluctuations at a given energy density decreases 
with the number of particles. In other words, given a 
number of different initial equilibrium conditions, only a 
fraction n of them will give rise to breathers. To give 
a quantitative idea, for the FPU model with eo = 1, n 
is only about 50 % for a chain of N = 20 particles but 
rapidly approaches 100 % already for N > 100. In the 
absence of spontaneous localisation, the non-linear sys- 
tems show no difference with respect to the linear ones 
and E(t) decays according to Eqs. . 



C. The effective— exponent analysis 

We have repeatedly mentioned that for the FPU model 
no stretched exponential relaxation is observed, even in 
presence of a weak on-site force. In this section we wish 
to remark this statement by reporting a more detailed 
analysis of the energy decay. This can be accomplished 
by defining an effective exponent through the loga- 
rithmic derivative 



C(*) = 



d[\nV(t)} 
d[\nt] 



(7) 



It is clear from definition (^) that if a portion of the en- 
ergy decay curve goes as E{t) = E(0) exp[— (i/ r )°1 this 
would result in a plateau of height a in the correspond- 



ing £(t) curve. For the problem at hand, interpretation 
of numerical data with this indicator is however a deli- 
cate matter. Actually, we want to remark that a naive 
analysis can lead to misleading or even incorrect conclu- 
sions 

The most convincing way to illustrate this is to discuss 
the case of the harmonic chain where localisation and 
stretched-exponential behaviour are obviously excluded 
a priori. From Eq. (Q), it is easily seen that C should start 
from 1 and monotonically vanish as 1/ lnt for t 3> To. On 
the other hand, we already learnt that in a finite chain 
a further crossover to an exponential law exp(— 2t/r JV _ 1 ) 
occurs so that £ must approach again 1 for t 3> t n _ x . 
The net results of those competing trends is that, for 
a finite N, £ displays a broad minimum Cmin> at t ~ 
t jv „ 1 . The value £min (which is independent of 7) can 
be estimated by noting that from Eqs. (Q) and (Q) one 
has k l/ln[27rt/ro] for t ^ tq, and hence £ m i„ w 
l/21n[4A r /7r]. A simple plot of the curves shows how the 
the minimum Q m i n can be incorrectly interpreted as a 
plateau if observed on a too short time scale. Of course, 
the same scenario is likely to appear in relatively short 
FPU chains whenever the chosen initial conditions do not 
yield localisation (see Fig. ||). 

When localisation does occur, the residual state decays 
exponentially |(|. Therefore one expects C to converge 
to 1 as the residual state is approached. Actually, one 
observes first very small values of £, while such a con- 
vergence does occur only on much longer time scales. To 
see why, let us consider a one-breather residual state. 
The total energy decays as E(t) — Ebe~~ t / Tb , where Eb is 
the initial energy of the breather and Tf, its characteris- 
tic decay time. The time constant t& is huge: in fact, it 
is roughly inversely proportional to the breather ampli- 
tude at the edge sites of the chain. From definitions (0) 
and (§) one has 



c 



t/n 



t/n - HE b /E{0)} 



(8) 



Since — ln[Eb/ E(0)} is a number of order 1, we see that, 
at times such that t <C tj, ln[J5&/J5(0)], C is very small. 
The numerical results in this case are just a deceptive 
byproduct of the wrong normalization. The exponent £ 
will start converging to 1 only when t ^> Tb\n[Eb/ E(0)], 
which is a huge time, being also E b <C E(0). In this case 
the C analysis may be misleading in identifying the nature 
of the decay law. The way around such problem would 
be to calculate the effective exponent by normalizing the 
total energy to Eb . However, this solution is in turn 
complicated by the intrinsic uncertainty in locating the 
time origin for the decay of the breather state. 



III. 2D LATTICES 

Since a careful study of the 2D Klein-Gordon model 
has already been presented [||, we focus here on the 2D 
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FIG. 3: Plot of for a short FPU chain (N — 20) for different 7s and eo = 2.. (a) Short linear time scale, (b) Long 
logarithmic time scale. The dotted line is the 7-independent minimum („,;„ (see text). Localisation does not even occur for 
this initial condition but one can incorrectly conclude in favour of a stretched-exponential law. 



FPU model. In this respect, the most natural extension 
of model (Q) should involve a two-component displace- 
ment vector. For the sake of simplicity, we rather con- 
sider here a scalar model with only one degree of freedom 
per lattice site. Actually, a preliminary series of sim- 



an N x N lattice, and therefore factorize, yielding 



ulations of the 2D vector version showed no appreciable 
differences with what reported here, and will be discussed 
elsewhere. 



With the same choice of non-dimensional units in- 
troduced in section ||, the model on a N x N lattice 
(i,j = 0,..., N— 1) with damping on all edges is defined 
by the equations of motion 



N-X 



V'(uij+i - Uij) - V'{uij - u itj -i) - r i'j Wj>,«) 

p,q=0 

where = j[gi, p 6 jtq + S itP g jt g - gi, p gj, q ], with g i>p = 
3i.p[3p,o + $p,N—i\- Since localisation has been shown to 
be strongly inhibited by fixed-ends boundary conditions, 
we consider here free-ends BC. 



E(t) 
£7(0) 



,-t/ro 



"7. 



-2t/r Q 



for t <C t 
for t > t 



27t(Vt , 

(9) 

with t q — N/2j. In particular, the law (H) applies in 
the absence of localisation, as well as during the first 
stage before the onset of localisation. In Fig. | we plot 
the energy decay curves of an FPU lattice with TV = 32 
relaxing from different values of the energy density eo- 
The first stage of the decay indeed follows the theoreti- 
cal prediction Eq. (||). Moreover, it is clear that locali- 
sation does not occur below a finite value of eo (located 
between eo = 0.2 and eo = 0.3 in this case). In par- 
ticular, in the absence of localisation the energy decay 
is again described by the two-crossover scenario. Notice 
that, also in this case, failure to spontaneously localise 
energy results in a minimum of the effective exponent £, 
due to the competition of the power law decay 1/t and 
the exponential law exp(— 2i/Vjv-i). 



Fluctuation— activated localisation 



Calculations identical to the one described in sec- 



tion II A can be easily extended to the case of a sim- 
ple 2D harmonic lattice. It turns out that the contribu- 
tions from the two spatial coordinates are identical for 



The relaxation process described in terms of the en- 
ergy decay curves is identical to the two-crossover pic- 
ture which applies to the ID FPU case. Nonetheless, 
the scenario markedly differs from the ID case in two 
respects. First of all, the mobility of the localised excita- 
tions which spontaneously emerge is drastically reduced. 
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FIG. 4: 2D FPU lattice, N = 32, 7 = 0.1. Plot of 2?(t) for 
different values of the initial energy eo (symbols) . The dashed 
line is a plot of Eq. (^) for t <C r . The inset shows the cor- 
responding effective exponents £(t) calculated from Eq. (Q). 
Notice the minimum associated with the absence of localisa- 
tion. 



In fact, after a first stage in which strong interaction is 
observed, they arrange on a "random lattice" , which, on 
the time scale of our typical simulations, appears indeed 
to be frozen. One of such states is illustrated in Fig. ^ 
(a). This is presumably due to a smaller "scattering sec- 
tion" in 2D. Obviously, due to the presence of dissipation, 
this state decays on a much longer time scale, which (as 
discussed above) scales exponentially with the linear size 
of the lattice. 

The second important difference with the ID case 
stems from the appearance of a finite energy threshold A 
for the existence of breathers (akin to envelope solitons in 
the small-amplitude limit) in 2D [jl4|. As a consequence, 
we expect that localised modes are generated in the relax- 
ation dynamics only by fluctuations that are large enough 
to overcome such threshold. The spontaneous excitation 
of breathers can thus be seen as an activated process. 
Accordingly, their number will be exponentially small in 
the ratio between A and some quantity measuring the 
strength of fluctuations. Therefore, one expects the av- 
erage density of breathers (i.e. the average number of 
breathers per lattice site) in the residual state to follow 
an Arrhcnius law of the form 



In our numerical simulations we clearly observe that 
n B strongly depends on the initial specific energy. In 
Fig. U (b) we plot (n B ) for two system sizes as a func- 
tion of eo, along with a fit performed with the law 
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(n B ) oc exp(-/3A) . (10) 

From Eq. (|lo|), it is tempting to identify (3 with some 
inverse temperature. However, one has to keep in mind 
that we are dealing with a non-equilibrium process and 
this identification can only make sense if energy release 
is adiabatically slow. If this is true, 1//3 should be pro- 
portional to (and smaller than) the initial temperature, 
which, in turn, is roughly proportional to the initial en- 
ergy density. 



(b) 

FIG. 5: 2D FPU lattice, N = 80, 7 = 0.1. (a) e = 1, sym- 
metrised site energies in the residual state state, (b) Average 
breather density (n B ) vs initial energy density for N = 30 
and N — 50 and Arrhenius plot. The inset shows the average 
density measured in the N — 50 system vs 1/eo in lin-log 
scale, and an exponential fit. 
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(n B ) — C exp(— A'/eo), with A' oc A. As it shows, the 
agreement is good. We conclude that the spontaneous lo- 
calisation of energy in our system is indeed an activated 
process. In particular, the fit gives A' sw 0.9. We note 
that a quantitative analysis of the time to required to 
reach the residual state (localisation time) reveals that 
to oc l/(n B ) oc exp(/3A) [ p~5| | . By fitting the numerical 
data for the localisation time, it is possible to get an 
independent confirmation of the thermal activation sce- 
nario. In particular, we get A' 0.7, in good agreement 
with the value obtained above. 



B. Statistics of breathers' energies 

It is important to realize that the "pseudo-stationary" 
distribution of breathers in the residual state emerges as 
a result of two competing mechanisms, namely the birth 
out of fluctuations from the homogeneous state and the 
decay due to both the coupling with the external reser- 
voir and to the various inelastic interactions (phonon- 
breather, breather-breather etc). Although the former 
can be safely neglected, it is very difficult to understand 
the latter in full detail. Indeed, to the best of our knowl- 
edge, only a few studies are available |lf| [l7| and a 
full description in terms of elementary interactions seems 
hardly feasible. Nonetheless, simple statistical arguments 
can be of help in understanding the energy distribution to 
a greater detail. Let us consider the fraction of breathers 
■p(e, t) having energy between e and e + de with e > A (in 
the FPU model breathers of arbitrarily large energy ex- 
ist due to the unboundedness of the interaction potential) 
and let us define B(e) and D(e) to be some phenomeno- 
logical birth and decay rates, respectively. Furthermore, 
let us assume that birth of a breather can only occur by 
a spontaneous fluctuation from the ground state (e = 0) . 
Stationarity requires the flux-balance condition 



B(e)V(0) = D(e)V(e) 



(11) 



Since the birth is an activated process we expect B(e) oc 
exp[— 0e] pl|. This assumption yields (for a constant 
V(0)) 



no) 

D(e) 



exp(-/3e) 



(12) 



In this simplified description, the prefactor of the ex- 
ponential term is thus interpreted as a measure of the 
effective breather lifetime r(e) (£>(e) oc l/r(e)) and is 
a priori unknown. In all simulations we observe that 
small-amplitude breathers decay more easily. This ob- 
servation agrees with the accepted existence of a pre- 
ferred energy flow from small-amplitude breathers to 
large-amplitude ones in breather-breather collisions jl7) . 
Moreover, Eq. (|l|) requires that t(A) = 0. Thus, we 
postulate 



D(e) oc (e-A)" 



(13) 



where the exponent z > can be estimated by measuring 
the distribution of breather energies in the vicinity of 
the threshold A. We can estimate the average den sity of 
breathers in the residual state from Eqs. (|l2] ) and (^3|) as 



) de oc (3 



-(z+l) e -f3A^ 



(14) 



This result is consistent with our initial hypothesis 
(Eq. ©). 



From Eqs. (|12|) and (g_3|) it follows that the average 
breather energy (e) can be expressed as a function of (3 
and A. Hence, we can write the normalised distribution 
Eq. ( |l2| ) as a function of the three parameters z, (e) and 
A in the following fashion 



V(e) = 



(* + l) 



z + l 



r(* + l)«e)-A) 



(e)-A 



A 



(15) 

where T(z) is the gamma function and (e) = A+(z+l)//3. 
From the point of view of the numerics, it is more accu- 
rate to deal with the cumulative (integrated) distribu- 
tions. These are easily obtained from the histograms 
without actually performing the integration, by not- 
ing that they are nothing but rank-size plots with the 
axes inverted. The cumulative distribution of the func- 
tion (|l|) is 



C(e) 



V(e')de' 



1 



r(z + i) 

x 7 (z + l, (z + l) 



A 



(e) 



(16) 



where 7(2, e) is the incomplete gamma function, defined 
as 



7 (z,e) = / y z L e v dy 



(17) 



In principle one could evaluate A explicitly by calcu- 
lating the energy of exact breather solutions of vanish- 
ing amplitude |Q. However, we expect that A <C (e). 
Therefore, if the population of breathers used to calcu- 
late the experimental distribution C(e) is large enough, 
we can assume that the smallest breather energy recorded 
£min is close to the threshold A. Hence, we can per- 
form a two-parameter fit with the theoretical prediction 
Eq. ( p^| ) with A = e m i n . The quality of such fits is il- 
lustrated in Fig. ^ for a lattice with N = 70. As it 
shows, the agreement of our simple model with the nu- 
merics is excellent. In particular, the analysis of the 
small-energy portion of the distributions yields z = 2 
(see inset in Fig. ^). We stress that the result of our 
fits were well reproduced irrespective of the particular 
realisation of e min . The fit of the experimental distri- 
butions over the whole range of available energies yields 
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The analysis of energy distributions also provides 
an independent confirmation that localisation is in the 
present case an activated process. To check this, we fit- 
ted the distributions obtained by letting systems of given 
size relax from different values of eo- Then, we evaluated 
P from the best-fit values of the fitting parameters z and 
(e) according to the relation 



1 





z + l 



(18) 



As shown in Fig. ^, j3 turns out to be inversely propor- 
tional to eo as expected. 



IV. CONCLUSIONS 



FIG. 6: 2D FPU lattice, N = 70, 7 = 0.1 and e = 1. Exper- 
imental cumulative distribution of breather energies (dashed 
line) and fit with formula ([Ti]) (solid line). Both curves have 
been shifted to the left by e m i„. The breather population is 
here of 1742 breather "events", recorded in 50 different real- 
isations of the initial condition. The inset shows a quadratic 
fit of the small-energy portion of the distribution. 



a slightly greater value, z w 2.6. This is because the 
large-energy tail of the distributions accounts for the 
rarest events, which are under-represented within rea- 
sonably large breather populations [jl9| . In order to con- 
vince oneself that this is the case, it is enough to in- 
troduce an energy-dependent weight function in the fit, 
which gauges the relative weights of the small-energy and 
large-energy portions of the curves. In this case, the best 
estimates of z monotonically decrease towards the value 
z = 2 upon lowering the relative weight assigned to the 
large-energy region. 



80 - 
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40 - 



20 - 




2 4 6 8 10 



FIG. 7: 2D FPU lattice, 7 = 0.1, N = 50. Plot of Eq. @ 
evaluated from the best-fit values of the floating parameters 
z and (e) vs initial energy density (symbols) and linear fit. 



The mechanisms yielding the spontaneous formation of 
localised periodic excitations, i.e. breathers, in spatially 
discrete non-linear dynamical systems are of primary im- 
portance for the understanding of their physical interest. 
In this paper we have focused on the phenomenon of re- 
laxation to breather states by cooling from the bound- 
aries ID and 2D lattices. We have provided a detailed 
description of the many facets of this phenomenon by 
combining numerical studies with theoretical arguments. 
In particular, we have shown that in models with suffi- 
ciently weak "substrate" forces the energy loss process 
can be explained on the basis of a simple perturbative 
analysis, that applies independently of the boundary con- 
ditions and of the nature of the nonlinearity. In this 
sense, we can claim that this is a general feature of this 
class of models, where an initial exponential decay is fol- 
lowed by a power-law behaviour. In the second stage 
of the decay process, boundary conditions play a crucial 
role in allowing for the formation of a long-living breather 
state, when they are taken free. On the contrary, fixed 
ends yield complete cooling of the lattice eventually ruled 
by the exponential decay rate of the longest-lived Fourier 
mode. We have also pointed out that the presence of a 
sufficiently strong "substrate" force can turn the initial 
exponential decay to a stretched-exponential law, while 
maintaining all the other features of the previously de- 
scribed scenario. A theoretical explanation of the origin 
of this stretched exponential behaviour and its depen- 
dence on the strength and on the nature of the local 
potential is a problem that deserves further and more re- 
fined investigations. We have also commented about the 
technical difficulties that can be encountered for identi- 
fying the correct time behaviour in these models, where 
the crossover between different regimes can be easily mis- 
taken for an indication in favour of other scaling laws. 

Another crucial aspect that we have widely analysed 
concerns the main differences between ID and 2D sys- 
tems. At variance with the former case, in 2D we have 
found numerical evidence of the existence of a finite en- 
ergy threshold for breather formation. We have also 
proposed a statistical model that provides an effective 
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agreement with numerical results concerning the breather 
residual state. This state eventually sets in as an almost 
static stationary configuration, where breathers of differ- 
ent amplitudes coexist. The statistical model is based 
on the hypothesis that this state emerges as a result of 
the competition between breather activation and mutual 
interaction. Moreover, it allows to determine an effective 
breather lifetime and its dependence on energy. In this 
framework the initial energy density plays the role of the 
control parameter regulating the strength of fluctuations. 

In our opinion these results provide a satisfactory un- 
derstanding of the phenomenon of spontaneous breather 
formation by cooling. As a final remark, it seems worth 



investigating experimentally the capability of real sys- 
tems to store energy in the form of long-lived localised 
excitations. 
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